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Abstract 

In many applications, and in particular in seismology, realistic propagation media dis- 
perse and attenuate waves. This dissipative behavior can be taken into account by using 
a viscoacoustic propagation model, which incorporates a complex and frequency-dependent 
viscoacoustic modulus in the constitutive relation. The main difficulty then lies in finding an 
efficient way to discretize the constitutive equation as it becomes a convolution integral in 
the time domain. To overcome this difficulty the usual approach consists in approximating 
the viscoacoustic modulus by a low-order rational function of frequency. We use here such an 
approximation and show how it can be incorporated in the velocity-pressure formulation for 
viscoacoustic waves. This formulation is coupled with the fictitious domain method which 
permit us to model efficiently diffraction by objects of complicated geometry and with the 
Perfectly Matched Layer Model which allows us to model wave propagation in unbounded 
domains. The space discretization of the problem is based on a mixed finite element method 
and for the discretization in time a 2nd order centered finite difference scheme is employed. 
Several numerical examples illustrate the efficiency of the method. 

1 Introduction 

Real media attenuate and disperse propagating waves |19j . Our aim in this paper is to develop 
a numerical method to model such dissipative phenomena (dispersion plus attenuation) in the 
time domain. To do so we consider the linear viscoacoustic equation which is a convolution in the 
time domain, the viscoacoustic modulus being frequency dependent. Therefore, incorporating 
any arbitrary dissipation law in time-domain methods is in general computationally intense. The 
usual way to overcome this difficulty is to approximate the viscoacoustic modulus by a low-order 
rational function ^JJ El E2j This leads to replacing the convolution integral by a set of 
variables, usually referred to as memory variables, which satisfy simple differential equations 
that can be easily discretized in the time domain. 

Several methods have been proposed in the literature for incorporating realistic attenuation 
laws (e.g. frequency-independent or weakly frequency- dependent viscoacoustic modulus) into 
time-domain methods ^JIEHEHEIEUE]- We focus our attention in this paper on the methods 
proposed by Day and Minster (1984), Emmerich and Korn (1987), and Blanch, Robertson and 
Symes (1995). All three methods use some approximation of the viscoacoustic modulus by a low- 
order rational function. The first approach is based on the standard Pade approximation. The 
coefficients of the rational approximation are thus in principle known analytically. Numerical 
results obtained using this method show that the approximation is poor and the method provides 
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satisfactory results only for relatively short (in terms of the wavelength) propagation paths. The 
second approach is based on the rheological model of the generalized Maxwell body, which gives 
a physical meaning to the coefficients of the rational approximation. They are interpreted as 
the relaxation frequencies and weight factors of the classical Maxwell bodies, which form the 
generalized Maxwell body. This method provides good numerical results for long propagation 
paths, but some parameters, namely the relaxation frequencies are semi-empirically determined. 
Finally, the third method is based on the observation that for the frequency-independent case 
and for weakly attenuating materials the weight factors are only slowly varying and can be 
approximated by a constant. This method provides good numerical results, but also involves a 
semi-empirical choice of a parameter. 

Although, the previous methods give satisfactory results in the case of weakly-attenuating 
materials they fail in media with large attenuation. This case was considered in a recent paper 
PQ, where the authors propose an analytic method for computing the best (optimal) rational 
approximation for the frequency independent case. They also propose a generalization of the 
algorithm presented in ^H] which leads to very good results in the case of highly attenuating 
media and a frequency- dependent viscoacoustic modulus. 

After a brief overview of the basic theory describing wave propagation in viscoacoustic media 
(section EJ) , we describe in section 01 the approximations proposed in CHI and [TO] . 

Considering long propagation paths, we test the performance of the different approximations 
and find that the best method, using the smaller number of unknowns while providing satisfac- 
tory numerical results and involving the least number of empirically determined values, is the 
one proposed by Emmerich and Korn (1987). We thus chose this method for approximating 
the viscoacoustic modulus. Note that a slight variation of the method proposed in ^Hj is used 
here, based on a different way of distributing the relaxation frequencies in the bandwidth of the 
incident pulse. 

In section El we incorporate this approximation in the velocity-pressure formulation for vis- 
coacoustic waves. Our choice of using the first-order-in-time system of equations, instead of the 
more classical second-order one, is motivated by the use of the fictitious domain method and 
the perfectly matched absorbing layer technique. In P] the authors proposed a similar approach 
using the mixed velocity-stress formulation for modeling wave propagation in viscoelastic media. 

The fictitious domain method (also called the domain embedding method) has been devel- 
oped for solving problems involving complex geometries [21 E21 EH1 , and, in particular, 
for wave propagation problems ^Hl EDI EH E] ■ In the framework of seismic wave propagation 
we apply this method to model the boundary condition on the surface of the earth (section EJ). 
Its main feature is extending the solution to a domain with simple shape, independent of the 
complex geometry, and to impose the boundary conditions with the introduction of a Lagrange 
multiplier. Thus, the solution is determined by two types of unknowns, the extended unknowns, 
defined in the enlarged simple shape domain and the auxiliary variable, supported on the bound- 
ary of complex geometry. The main advantage is that the mesh for computing the extended 
functions can now be chosen independently of the geometry of the boundary. 

The Perfectly Matched Layers (PML) technique was introduced by Berenger [HI E] for 
Maxwell's equations and is now the most widely-used method for the simulation of electro- 
magnetic waves in unbounded domains (cf. (HU03EE]) It has also been extended to the case 
of anisotropic acoustic waves isotropic [2H1 and anisotropic elastic waves ^ElE]- This tech- 
nique consists in designing an absorbing layer, called a perfectly matched layer (PML), that has 
the property of generating no reflection at the interface between the free medium and the arti- 
ficial absorbing medium. This property allows the use of a very high damping parameter inside 
the layer, and consequently of a small layer width, while achieving a near-perfect absorption of 
the waves. We apply here the PML model in the case of viscoacoustic waves (section EJ) . 
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Another advantage of the first-order formulation over the second order one, is that it is easier 
to implement in heterogeneous media, since it does not require an approximation of spatial 
derivatives of the physical parameters. To discretize this formulation in space we use a mixed 
finite- element method which is a modification of the method proposed in pjj. More precisely, in 
[H] the authors designed new mixed finite elements, the so-called Qjfe — Q k elements, inspired 
by Nedelec's second family which are compatible with mass lumping, and therefore allow 
to construct an explicit scheme in time. A non-standard convergence analysis of the Qf+i — Q k 
elements was carried out in 0. However, numerical results obtained recently (cf. show that, 
when coupled with the fictitious domain method, these elements do not provide satisfactory 
results. This is why we use here instead the (jf+i — P k+1 elements for which convergence of the 
fictitious domain method was obtained 

To show the efficiency and robustness of the method we present in section several numerical 
results. In particular, numerical and analytical results are compared and good agreement is 
obtained between the two. 



2 Viscoacoustic wave propagation 

In an isotropic viscoacoustic medium occupying a domain Q G M. d , d = 1,2,3, the relation 
between the pressure p(oS) = p(x, uS) and the displacement u(u>) = u(x, u) in the frequency 
domain is, 

p(oj) = fj,(uj)divu(u). (1) 

Here, /j,(uj) is the complex, frequency-dependent, viscoacoustic modulus. 

The dissipative aspect of a material is often described by the quality factor Q, defined as the 
ratio of the real and imaginary parts of the viscoacoustic modulus. It expresses how attenuating 
a material is and corresponds to the number of wavelengths a wave can propagate through the 
medium before its amplitude has decreased by e _7r , 

where <p{u) is the phase of /j,(uj). 

In seismic applications, Q is usually assumed to be frequency- independent or only slowly 
frequency-dependent. In this case (i.e. when Q is constant in frequency), the viscoacoustic 
modulus is given analytically by Kjartansson's model |2(i| . 

/ fa \iarctaniQ- 1 ) 
A»M = Mre/ ■ (3) 

This analytical formulation will be useful for validation of the numerical results in the next 
sections. 

In the time domain, the constitutive relation £Q) is expressed in terms of a convolution 
operator, denoted here by *t, 

p(t) = n(t) * t divu(t). (4) 

The discretization of this equation requires saving in memory the whole history of the solution 
at all points of the computational domain and is thus very expensive. To overcome this incon- 
venience, we approximate the viscoacoustic modulus by a rational function in frequency, as was 
proposed in ^3 El H21 Ei]- It is convenient in the following to introduce the relaxation function 
R(x,t), defined by, (see Figure 

/x(x, t) = ; i?(x, t) = L R (x) + Mx) J o + °° r(x, u/K^d/) H(t), (5) 
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where hr is the relaxed modulus, 



W?( x ) = lini R (t)> 



Hu is the unrelaxed modulus, 



/i[/(x) = (ir(x) + fyt(x) = lim J?(i), 

r(x, u/) is the normalized relaxation spectrum satisfying 

y+oo 

/ r(x, uo')dw' = 1, 
•/ 

and is the Heaviside function. 

R(t) 

Su 



t=0 



Figure 1: Schematic example of the relaxation function R(t) 
Using /i(x, t) defined by (0) in Q gives, 

/t y+oo 
I u/r(x, u')e~^ ( t ~ T ^div\i(-x.,T)du)'d,T. 
-oo ^0 

We now assume that the relaxation spectrum can be discretized by L single peaks of ampli- 
tude ai at relaxation frequencies cjj, I E [!•••£], 



( x >^) = ^2oh(x)S(uj - Wj(x)) ; ^az(x) = 1, 



z=i 



In this case, we get, 



i?(x,i) ~ R t (x,t) = ^(x) +^a,(x)e- w 'W*j 



and 



M(X,W) ~ W(X,W) = /Zfl(x) 1 + > ; ■ — . 

In ©, we introduced yz(x) defined by, 

y/(x) = — — a;(x), with the normalization relation yi (x) 



(6) 



m x 



M x ) 
w?( x ) 



Notice that equation © can be obtained if one assumes that /x(x, to) can be approximated by a 
rational function of (iw), 

fj,{x,u) = , v (7) 

with Pl and Ql being polynomials of degree L in [ioS). Then © can be interpreted as an 
expansion of Q into partial fractions [^j. Thus approximating the viscoacoustic modulus by a 
rational function is equivalent to approximating the relaxation spectrum by a discrete one. 

For computational reasons, it is natural to search for rational function approximations of the 
viscoacoustic modulus, which minimize the ratio: number of unknowns/accuracy. We therefore 
address in the following the question of finding an accurate low-order approximation of the 
viscoacoustic modulus. 

3 Approximation of the viscoacoustic modulus 

We now briefly introduce the different approximation methods previously proposed in the 
literature. 

3.1 Pade approximation method 

The use of the simple Pade approximation in the framework of viscoacoustic wave propaga- 
tion was proposed in WJ\. Letting z = — — and introducing, 

X(x, z) = / — —du> , 

Jo 1 - u'z 

h(x,oj) can be re- written in the following form, 

/i(x, lo) = /U[/(x) + <fyi(x)zx(x, z). 

The Pade approximation is then used for expanding x(x, z) into a rational function with nu- 
merator of degree L — 1 and denominator of degree L. Using the well-known (jSI3^B) relations 
between Pade approximations and orthogonal polynomials one gets, 



(v) = Errr' 

^-^ 1 — uji(x)z 



l=1 - 

where u>i are the zeros of the orthogonal polynomial Pl, and A; are the residuals given by, 



kL-iPL-M^p^i^y 



kL being the leading coefficient of Pl and where the prime denotes the derivative of Pl- Recall 
that the orthogonal polynomials are defined by, 

P n (uj')P m (uj')uj'r(u;')duj' = 5 mn , 

m 

where 6 mn is the Kronecker symbol. When the quality factor is constant over a frequency 
band, A; and uji can be obtained in closed form. Moreover, when Q < 1, the relaxation spectrum 
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r(x, u}) is proportional to lu , Assuming that r(x, oj) is zero outside the frequency interval 
[fix , Q2] we obtain the approximation, 

f, n2 ~ m ^( x ) \ 

where u>i = -[xi(Q2 — 01) + 172 + 0,1], x\ and v\ being respectively the zeros and weights of 
the Legendre polynomials. Notice that the relaxation frequencies are in this case equidistant 
on a linear scale. The main advantage of this approximation is that all data are analytically 
determined. For more details on this method the reader can refer to |17j . 

3.2 Generalized Maxwell Body approximation method 

We describe here the method proposed in ^H]- First let us re- write ijHJl as, 

W (x, u) = /i R (x) + frt(x) V . "^ X) 7 . (9) 

Each term of (El) can be interpreted as a classical Maxwell body with viscosity a\ — and elastic 

modulus ctiSfi. The term hr in (JHJ) represents an additional elastic element. The Q-law for the 
generalized Maxwell body approximation can be obtained from (J§J), 



EL / \ u>i (x) 



Assuming now that <5// -C /ijj, fflTTfl becomes, 



L 



This means that Q(w)^ 1 is approximately the sum of n Debye functions with maxima a>i 

2//R 

located at frequencies u>i. If Q is fairly constant in a frequency band, the most natural choice for 
the relaxation frequencies u>i is a logarithmic equidistant distribution. In this case, to obtain a 
good approximation of the distance between two adjacent relaxation frequencies should 

be chosen smaller or equal to the half- width of the Debye function (1.144 decades). In [T8] 
two ways for choosing u>i were proposed: ui can be chosen logarithmically-equidistant in the 

frequency band [01,02] or determined by loi = — d °, m where uidom is the dominant (central) 

10' 

frequency of the source considered in the simulations. In both cases, the coefficients m are 
obtained by solving the overdetermined linear system 

A ^(x) - Q- 1 (x,5 fc (x))5 fe (x) ~ x _ 

> I/j(x)w fc (x) , .„ „ , .„ = Q (x,o;fe(x)), A e 1,2, ..AT , 12 

i=l 

where, uik are defined by 



c^(x) 2 +cj fe (x) 2 
wi = Oi 



,02. 



1 



Let us remark that the determination of u>i for this approximation is based on an empirical 
study. 
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3.3 The r-method 



This method, proposed in [1Q] is based on the observation that dissipation due to only one 
"Maxwell Body" can be determined by a unique dimensionless parameter r. More precisely, for 
Q > 1 and L = 1, equation ifTTjl becomes, 



Q{x,u) - 



T X 



l + ( 



wi(x) ■ 



where r = y\ <C 1. It is then easy to see (cf. |10|). that wi essentially determines the 
frequency behavior of Q while r determines its magnitude. In the general case for L > 1, and 
when one seeks an approximation of a constant Q value, yi are quasi- const ant 

and equation 111 111 can be approximated by, 



, =1 i + 



In (fTTTf) . Q(o;) 1 is linear in r. One can therefore find the best approximation, in the least-squares 
sense, over a predefined frequency range to any Qq by minimizing over r the expression, 

J= / {Q-\u,u h r)-Q^) 2 dw. (14) 
./fii 

The approximation of the viscoacoustic modulus in this case is, 

«^ = "«w( 1 + Sl7^S5o)- (15) 

The relaxation frequencies u\ are chosen, as for the "Generalized Maxwell Body" method, equidis- 
tant on a logarithmic scale. Equation (|TK|) leads in general to an over-estimation of the value 
of Q. Thus the authors in ^2] suggest to use in the definition of J ifTl]) a value for Qo slightly 
smaller than the desired one. This value is also chosen empirically. 

3.4 Comparison of the different approximation methods 

To test the accuracy of the different approximation methods previously presented, we com- 
pute the response of a one-dimensional viscoacoustic homogeneous medium to the following 
pulse, 

( 2rrt\ /47ri\ 
s(t) = sin ( — J - 0.5sin ( — J for < t < T, T = 0.3s. (16) 

The solution is obtained by convolving the source function s(t) with the dissipation operator 
D(t) (the Green's function for the ID problem). For an arbitrary dissipation law, the Fourier 
transform D(uS) of D(t) is given by ^ 



D(,) = e" i * ) ( 1 -^), (17) 

where c(u> r ) is the phase velocity at the reference frequency u> r , v{oj) the complex velocity, and 
t* = — — r— -. — r the dissipation time. For a frequency independent Q, the value of — j^— = 

C{UJ r )Q{LO r ) ' * U \ UJ ) 

H^r)| 
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Figure 2: Comparison between the different approximation methods. Correlation coefficient 
between the exact and the approximated solutions as a function of the dissipation time. 



can be determined from equation (0J) combined with one of JEJ or (|lf>p. depending on 
the approximation method used. 

In the numerical example, we want to approximate Q = 20 over the frequency range 
[10" 2 ,10 2 ] Hz, likein [IH]. 

To better illustrate the results, we present in Figure El the evolution of the correlation coef- 
ficient between the exact solution (the one obtained for the viscosity modulus calculated from 
Q) and the different approximated ones (calculated with the viscosity modulus provided by 
JSj, JEJ or ifTKjl ) as a function of the dissipation time. 

More precisely in Figure we compare the results obtained with the following approxima- 
tions, 

• Pade approximation with L = 5. 

• Maxwell Body approximation with L = 3 and relaxation frequencies chosen logarithmically- 
equidistant over the frequency range [10~ , 10 1 " 5 ] Hz (cf. (IHj)- We call this choice 
method 1. 

• Maxwell Body approximation with L = 3 and relaxation frequencies chosen equidistant 
on a logarithmic scale, such that, oji = 2 ^° m (cf. [IB])- We call this choice method 2. 

• The r-method with Qo = 17.6 (value proposed in ^0] to model the propagation in a 
viscoacoustic medium with Q = 20) and L = 3. 

The results illustrated in Figure show that the Pade approximant provides good accuracy 
only for short dissipation times, as demonstarted in [18J. The r-method provides a good ac- 
curacy/number of calculations ratio. However, we did not choose this method because Qo has 
to be calibrated empirically in order to get good results. The "Generalized Maxwell Body" ap- 
proximation method seems to be a good compromise between accuracy, number of calculations, 
and implementation simplicity. As our aim is to simulate viscoacoustic wave propagation in 
heterogeneous media for large dissipation times, we chose a method which is a hybrid of the 
Maxwell approximation methods 1 and 2. 
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3.5 Proposed method 



In practice, the source type used depends on the application of interest. In our case, the 
main application of interest is seismic wave propagation for which a Ricker wavelet is often used 
as source function, 

f(t) = -2a 2 ^1 - 2a 2 - yj ^ exp (^-a 2 - y^j , for < t < y, with a = vr/ . (18) 
In Figures [HI and @] we display the source function l[T8|l and its spectrum for two different 




0.5 1 1.5 2 

t(s) 



Figure 3: The Ricker wavelet f(t) for /q = 
2.5Hz (solid line) and for fo = 1Hz (dashed 
line). 




v(Hz) 

Figure 4: The frequency spectrum of the 
Ricker wavelet for fo = 2.5 Hz (solid line) and 
for /o = 1Hz (dasged line). 



values of fo. Compared to the source function defined by lfTB|) . the Ricker wavelet has a broader 
frequency spectrum and method 2 did not give as good results in this case as the ones obtained 
with the source (|T6|l . Following the ideas in ^H], we want to find a way to choose the frequency 
band [01,J72] as a function of the source type and then determine the relaxation frequencies u\ 
logarithmically equidistant in this bandwidth. We found that a good choice for a Ricker wavelet 
type of source is [Ql, Q2] = [ , ui max ] , where u> max is the maximal frequency of the employed 

source spectrum. None of the above approximation methods is completely satisfactory in our 
opinion because the choice of the relaxation frequencies is always empirical. 

To avoid this, one can follow the approach proposed in £Q where a non-linear minimization 
problem is considered which permits to determine all the coefficients (both yi and u>i V/ G [1, L\). 
However, this method is more expensive and although it improves the accuracy of the solution 
for media with high damping (Q < 10) it provides quite similar results with the proposed method 
for propagation in weakly attenuating media (Q > 10) ^p. As we are interested in media with 
quality factors greater than 10, we will use in the following the linear minimization method 
(system (|T2"|l h 
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4 The mixed velocity-pressure formulation 

By incorporating (j^J into Q we get, 

p(x, w) = /i R (x)div(u(x, w)) + /i/?(x) V] — rT div(u(x, u)). (19) 

^ lW + UJl{X) 

We now introduce the memory variables rji denned by, 

(iu + ui(x))rn(x,u;) = /x R (x)yz(x)div(v(x, u)), (20) 

where v is the velocity, i.e., the time derivative of the displacement u. Equation (fflTj) in the time 
domain becomes, 

ftftOM) + m( x ),K(x,t) = ^ fl (x)y,(x)div(v(x,t)). (21) 
Using the definition of 77/ and multiplying lfT?I)l by (iu;), we get, 



(iw)p(x,w) = j u jR (x)div(v(x,u;)) + ^(iu;)r#(x 
or equivalently in the time domain, 



!=1 



t = — > + £f- (22 , 

Combining (|22|l. ff2Tf) and the equation of motion, we obtain our final system of equations, 

( dv 

P^"Vp = f inOx]0,r], 
i=l 

^ +ujirji = fi R yidiv(v),yi in Ox]0,T]. 

Equivalently, one can chose to eliminate the pressure and obtain a second-order-in-time equation 
for the displacement by introducing adequate memory variables |18j . We prefer, however, the 
first-order velocity-pressure formulation for the following reasons, 

• It can be coupled with the fictitious domain method for taking into account diffraction by 
objects of complicated geometry. 

• A perfectly matched layer model (PML) can be written for this system. This permits us 
to simulate efficiently wave propagation in unbounded domains. 

• This system is easier to implement in heterogeneous media, since it does not require an 
approximation of the spatial derivatives of the physical parameters. 

An equivalent first-order velocity-pressure system is proposed in ^2] and ^0]- In ^2] the authors 
used a pseudospectral method for the discretization while in a staggered finite difference 
scheme was used. Our aim being to couple this system with the fictitious domain method, 
we propose here instead the use of a mixed-finite element method on regular grids. A similar 
approach was proposed in [3j where the authors use a mixed-finite element method to discretize 
the velocity-stress formulation for viscoelastic wave propagation. 
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5 Discretisation 



A mixed formulation associated to equations (|28|) is given by. 



' Find (v,p,iJ) :]0,T[i — ► X x M x (M) L s.t. 
d 



dt 



(pv,w) +6(w,p) = (f,w), 



Vw G X, 



0, G M, 



d PH ^ Mil 

4(— »7i,9) + (—rii,q) -6(v, 9) =0, V/, G M, 
t di ^/jy; 



(24) 



where H is the L-dimensional vector with components rji, and 



&( w ) Q) = Q divw dx, V(w, g) G X x M. 



The functional spaces are X = H(div; Q), and M = L 2 (Q). 

We now introduce some finite element spaces X^ C X, and C M of dimensions iVi and N2 
respectively. The semi-discretization in the space of problem {21}) is, 



(V h ,P h ,H h ) G L 2 (0,T;IR N i) x L 2 (0,T;ffi^) x L 2 (0iT; (j^L) g t _ 
= ^h, 

dP /( ^„ 



M v ^ + B h P h 



M 



M 



dt 

d(H, 



Y, M v 



1=1 



dt 



h)l 



B T h V h 



+ M u {H h ) l -BlV h = ^ VZ 



0. 



(25) 



dt 



where denotes the transpose of B^. 

In practice, we only consider regular domains in IR d , d = 1,2 that can be discretized with a 
uniform mesh 7^ composed by segments or squares of size h, depending on the dimension of the 
problem. The finite element spaces we use are 



X h = {w h in X I VK G %, w h \ K G (Qi) d } , 



and 



M h = {q h eL 2 /VKe T h , q h \ K G P (K)} . 

This mixed finite element was introduced in (Sj and is illustrated in Figure El 
When coupled with the fictitious domain method, this choice of finite elements presents some 
inconveniences. In particular, for the acoustic wave equation problem we cannot prove the 
convergence of the method from the theoretical point of view. Moreover, numerical results show 
that the method converges under restrictive conditions on the discretization parameters. Thus, 
when the method is coupled with the fictitious domains, we replace Mh by M\ defined by, 



Ml 



{q h £L 2 /\/K€T h ,q hlK G P X {K)} . 



This finite element is presented Figure El In this case convergence for the acoustic waves problem 
coupled with the fictitious domain method was obtained 
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Figure 5: Finite element v h G X h , (p h , (r]i) h ) G M h x M h 



For computational reasons, however, it is natural to seek a discretization which uses the least 
number of variables. In the proposed scheme pressure is thus discretized on the space M\ and 
the memory variables rji are discretized on M/,. The semi-discretization of the problem l(2"3|) in 
this case is, 

f (V h ,P h ,H h )e L 2 (0,T;ffl Af >)xL 2 (0 1 T;ffi*)xL 2 (0,r;(ffl*)i) s .t. : 



M^ + BlP h = Fi 
M idPh_y M ad(H h )i 

P M ^ P 



Bl' T V h = 0, 



(26) 



dt ^ "~ p dt 
1=1 

- M °^sP + M » {Hh)l - B ^ Vh = °' v/ - 

In both cases (pressure discretized on Mh or M^), we use a second order centered finite difference 




1 < 











— - 


p 






■ 




■ 


1 »; 


p p 







Figure 6: Finite element v h G X h , (p h , (rji)h) G Ml X 
scheme for the discretization in time (here presented in the more general case with the pressure 



discretized on M^), 



yn+l 

Ml Vh 



At 

p«+| _ 
Ml h 



G IR N ' x IR N2 x (ffi^) L 

L 



p 



At 



( H h)i 



At 



B i h T vr i 



B 0,T K+1 



(27) 



o,vz. 



y At ' ~~ u 2 

The numerical scheme (|27|l becomes explicit in time when an adequate quadrature formula is 
used to approximate the matrix M*. Note that the other mass matrices (M°, Mp, My, and 
M°) are diagonal, since the spaces and M\ are composed of discontinuous functions. For 
more details on the quadrature formulas used we refer the reader to pj. 
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6 Stability and dispersion analysis 

For the continuous problem, the energy is defined by 



1 1 / L L \ L 1 

This quantity is positive (for y/ positive) and we have, 

That is, the energy decreases as a function of time, which expresses the dissipative nature of the 
problem. 

In the discrete case a stability analysis based on energy techniques permits us to show that 
the discrete scheme is stable under the following CFL condition (in homogeneous media and for 
both choices and M\ for the pressure discretization), 

^ m \B h t(x + Y.y\^ (30) 



4 p 



i=i 



rp 4 rp . 8 

with \\Bi B h \^—r in ID and \Bi B h \\ ^ —r in 2D. Note that these are the usual CFL conditions 

obtained in the non-dissipative case multiplied by ^1 + . 

Furthermore, the dispersion relation also presents a similar aspect. For the the continuous 
problem we have, 

For the discrete problem in ID we obtain, 

2 fuAA _ A 2 t c 2 ( . 2 fkA x \\ ( A 2iy ; tan(^' 



sin = — - — sin 



1 + V — { ~> (32) 
1 + ^A^ + 2itan(^)J' ^ 



and in 2D we get (for both choices Mh and M\ of the pressure discretization), 



„M _ ^ u,f^A +rf [S^>\] fi + E A 2i "T t ( tL ) • (33) 



sin 2 I — — I — — - — ! sin 



x ^ x A t ui + 2itan ^ 2 / j 

In figure El we have plotted the dispersion and attenuation curves as function of 1/N (N 
being the number of points per wavelength used in the discretization) for a plane incident wave, 
whose incident angle is or 7r/4 for the 2D case. Note, in particular, that the ID scheme is no 
longer exact as it is the case in non-dissipative media. Depending on the angle of incidence, the 
2D scheme may be more or less dispersive than the ID one. 

Demonstration and details of the calculations for the stability and the dispersion relations 
for the discrete problem are exposed in the Appendix 1X1 and IbI 
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Dispersion curves for 7r/4 



' 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 
1/N 

Attenuation curves for 7r/4 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 



Dispersion curves for 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 



Attenuation curves for 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 



Figure 7: Dispersion and attenuation curves for a plane incident wave, with incident angle or 
7r/4 for the 2D case. 
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7 The fictitious domain method 



To model the free-surface boundary condition on the surface of the earth we use the fictitious 
domain method which has been developed for solving problems involving complex geometries 
OESEHlEUEij, and in particular for wave propagation problems [T51 1231 1531 IH] . 

We follow here the approach proposed in (§]. Consider the viscoacoustic wave propagation 
problem in a domain with a complex geometry such as the one described in Figure El The initial 
problem is posed in O with the free-surface boundary condition, v • n = on T, 

dv 

p— Vp = f in f2, 



dt 

dp drji 

dt ^ dt 

i=i 

^ + uJirji = /iR2/jdiv(v), V/ in Q, 



^i?div(v) in Q, 



(34) 



v • n = 0, 

p = 0, 



on r, 

on T D . 




Figure 8: Geometry of the problem: on the left the initial domain f2 and on the right the 
extended domain C. 

The main idea of the fictitious domain method is to extend the solution to a domain with a 
simple shape, independent of the complex geometry of the boundary, and to impose the boundary 
condition in a weak way by introducing a Lagrange multiplier. Following this idea, we extend 
the solution (v,p,r]i) by zero in the domain C (which is here a rectangle, see Figure EJ- We 
denote (v,p,fji) the extended solution and have, 

[vn] r = => v G H(div, C), [p] T £ 0, [j»] r + 0. 

Thus, system (|8*l^l for the extended solution, can be written (in the distributional sense), 



9v 

'at 



Vp = f + [p]n<5r 



dp s-^ drji 

Tt ~ 2^ -at = ^ dlv(v) 
i=i 



in C, 
in C, 



dr)i 
dt 



+ uifji = p R yidiv(v), V/ in C, 



(35) 



v • n 



p = 0, 



on r, 

on T D . 
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In 1)3 5|) we have two types of unknowns, the extended unknowns, defined in the simple shape 
domain C and the auxiliary variable [p], defined on the boundary V. We introduce A = [p] 
as a new unknown defined on T. This unknown can be interpreted as a Lagrange multiplier 
associated with the boundary condition on T. The variational formulation of the problem can 
then be written as follows, 



' Find (v,p,H,\) :]0,T[i — > H(div;C) x L 2 (C) x (L 2 (C)) L x HV 2 {T) s.t. 
d 



dt 



(pv, w) + b(w,p) - 6r(A, w) = (f, w), 



d -(—p,q) %g) -b(v,q) = 0, 



dt n R 



T A—vi,q) + (—m,q)-b(v, q ) = o, v/, 

dt fiRVi n R yi 
. 6r(/i, v) = 0, 

where 



Vw G ff(div;C), 

VgeL 2 (C), 
e l 2 (C), 

G FV2(r), 



6r0i,w) = y fiwnds, V(/i,w)eff 1/2 (r)xff(div;C). 

For the discretization of this problem we consider a structured volume mesh 7^ on C, and 
an irregular surface mesh Qh a on Y. The main advantage of this formulation is that the mesh for 
computing the extended functions can now be regular while the surface mesh is irregular and 
permits a good and efficient approximation of the geometry (see Figure El- 
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Figure 9: The two meshes used in the fictitious domain method: a structured volume mesh on 
the domain C and an irregular surface mesh on T. 

To discretize the volume unknowns (v,p,H) we use the finite element method described in 
section El while for the Lagrange multiplier A we use piecewise linear continuous functions on 
Qh s , i.e., the approximation space is, 

Gh a = {nh s G # 1/2 (r)/ VS(segment) G G hs ,n hs \ s G P^S)} . 

To simplify the presentation, we considered in system |(HH the homogeneous Dirichlet bound- 
ary condition on the boundary Tp. When the domain is infinite we use the perfectly matched 
layer model which will be described in the following section. 
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8 The PML method 



The perfectly Matched Layer model was introduced by Berenger (HIE] for Maxwell's equations 
and is now the most widely-used method for simulating wave propagation in unbounded domains. 
The reader can refer to j^l^EH] for electromagnetic waves, to 0] for anisotropic acoustic waves 
and to ^ElE] for elastic waves. The popularity of this model is due to its simplicity and efficiency. 
Its most astonishing property is that for the continuous problem the reflection coefficient at the 
interface between the layer and the free medium is zero for all frequencies and angles of incidence. 
To derive the PML for the viscoacoustic system (|23|) we follow the approach proposed in ^2] 
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Figure 10: PML in the x-direction: the physical medium is on the left and the absorbing medium 
is a layer of width 5. 



which applies to any first-order linear hyperbolic system. We present here the construction of a 
PML in the x-direction (see Figure [Tf)|) . Deriving then the PML for the other boundaries and 
the corners of the computational domain is a straightforward application of the same technique. 

Following pH] we construct the PML model in two steps: (i) We split the solution (v,p, {rji}) 
into two parts (vll,pH, {w}) and (v ± ,p ± ,{r]j L }), with the parallel part being associated with 
the derivatives in the y-direction (direction parallel to the interface between the PML and the 
physical medium), and the orthogonal part associated with those in the x-direction. (ii) We 
introduce damping only on the orthogonal component of the solution. 

When applying the splitting step to 
we obtain, 



23|) by remarking that v" = (0,v y ) and v = (v x ,0), 



dp 
dy 



dv y 

dpW 
dt 

d Vi 
dt 



1=1 



dt 



dy 



(36) 



+ w/r/j 1 = fi R yi 



di 



dy 
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and 



dv x 
dt 



dp 1 - 

H 

dt 



dp 

dx 

L 



E djn __ uv x 
Fit ~ ^ R FIt 



dv x 



1=1 



dt 



dx 



(37) 



with 



(38) 



\ m = Tff +r]^, V/. 

To apply the damping on the orthogonal components it is simpler to consider system l(3*7jl in the 

d x (cf. 



frequency domain. Then the PML consists in replacing the ^-derivatives d x by 
[16]). Following this approach, system (|3*7|l in the frequency domain becomes, 

dp 



1uJ 



iuj + d(x) 



p{iu + d(x))dv x = — 
ox 



(iu + d(x))p- L - ^2(iuJ + d(x))r]l = ti R 



dv x 
dx 



(39) 



{in) (iu)(iu) + d(x))r]^~ + (iu + d(x))coirj^~ 



dv x 



where d(x) is the damping parameter which is equal to zero in the physical medium and non- 
negative in the absorbing medium. 

We now introduce new variables rji defined by, 



iuir]i = {iuj + d(x))r]^~, V7, 



(40) 



or equivalently in time domain, 



dm 

dt 



dm 



+ d(x)rrf; V/. 



dt ' ~ K ~ rH ' 

Using 1(30]) in l(39"|) and going in the time domain we get, 

dv x dp 
p— + P d{x)v x = - 



dp 1 

~dt 

drji 



+ d(x) P ± ~^Z^ = m 



1=1 
dv x 



dt 



dv x 
dx 



(41) 



dt 1 1 " ™" dx 

The final system of equations for the PML is iflTjl together with l(36|l . with p being defined by 
p = p" +p~ L . Note that the memory variables rji do not appear, and only the component rjj and 
the variables rji do appear, in this system. 

Using a plane wave analysis, it can be shown (cf. ^H]) that this model generates no reflection 
at the interface between the physical and the absorbing medium and that the wave decreases 
exponentially inside the layer. This property allows the use of a very high damping parameter 
inside the layer, and consequently of a small layer width, while achieving a near-perfect absorp- 
tion of the waves. Note that for a finite-length absorbing layer there is some reflection due to 
the outer boundary of the PML. 
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Remark 1 To discretize the PML we use the same scheme as for the interior domain. 



Remark 2 The damping d(x) is zero in the physical domain and non negative in the absorbing 
medium. In the numerical simulations it is defined as in [16], 



f _ 

I l°9 (7? V ' " 



for x < 
(I)™ forx>0 



(42) 



,RJ 25 

where R is the theoretical reflection coefficient, 5 the width of the PML and n = 4. 
In practice, we take R = 5.0 1CT 7 , and 5 » 30A X (depending on the wavelength). 



9 Numerical results 

9.1 Scattering from a circular cylinder 

In order to validate the proposed numerical method we consider in this section the canon- 
ical problem of a plane wave (Ricker wavelet) striking a viscoacoustic homogeneous circular 
cylinder. The geometry of the problem is displayed in Figure A homogeneous viscoacous- 
tic circular cylinder of radius a (domain O2) is surrounded by a homogeneous non-dissipative 
medium (domain Jli). We denote by Ti the interface between the two domains and ^2 ■ The 
physical characteristics of the media are p\ = 1000Kgr/m 3 , c\ = 3050m/s, Q\ = +00 in Q,\ 
and P2 = 1800Kgr/m 3 , C2 = 3050m/s and Q2 = 30 in f^- The source function used in this 
example is given by lfT%|) with fo = 2.5Hz. For this problem, the solution can be computed by 
an analytical method described in what follows. 

3500 m 



p =1000Kg.m " 
c,=3050 m.s 



n 2 

P 2 =1800Kg.m" 
Cj=3050 m.s"' 



Figure 11: The geometry of the problem: a homogeneous viscoacoustic cylinder of radius a 
(domain Q2) embedded in a non-dissipative homogeneous medium (domain Cli). 

Consider the following incident plane wave (with incident angle < 9 l < 2n), 

n=oo 

p\(x) = Ai ) e- in( - dl+ ^J n (k r)exp(ine); Vx = (r cos(8) , r sm(6)) G Q v (43) 



19 



Using the partial wave expansion we can express the solutions pj E Qj, j = 1, 2 in the following 
form, 

n=oo 

pi(x)=pi(x)+ ^ oin^ 1} (A;ir)exp(in0), VxGOi 

n=oo n= -°° (44) 

p 2 (x) = ai n J„(A;2(^)^)exp(in6'), Vx G Sl 2 , 

n=— oo 

with the first-kind Hankel function of order n , J n the Bessel function of order n and where 

— atan(-pi-) 



the wave number in f2 2 is computed by 

fc 2 (x,u;) = fe 2 (w) = 



U / 1U 



C-ref \^ref 



(45) 



To compute the coefficients a± n and &i n we introduce the expressions for p\ and p 2 , i.e., equa- 
tion in the transmission boundary conditions on Y\ (continuity of the pressure and the 
normal component of velocity). After projecting the resulting system onto the Fourier basis 

— exp(— imO); m G Z J we obtain, 



2vr 



din 



lljn{Xl)Jn{X2) ~ 72Jn(Xl)jn(X2) U -W^+f) 



7lM 1) (xi)^n(X2) + 72Jn(X2)^ ±J (xi) 



, HJn(Xl)H { n \xi) -llJn{Xl)H X n'{xi) A j -I»««+S) 

72«/n(X2)^ 1) (Xl)+7l^ 1) (Xl)«/n(X2) 



r(l)/ 



(46) 



with Z n (z) = dZ ^f^ , Xj = an d 7j = f 1 - The insertion of these expressions into ijHj) gives 
the final solution of the problem jHS|. Comparison of results between the analytical and the 
numerical solution are displayed in Figure El where we can see that good agreement is obtained 
between the two. 



Solution at point R\ 



Solution at point i? 2 



Solution at point R3 




Figure 12: Comparison between the analytical solution (dashed line) and the numerical solution 
(continuous line) at different observation points R\, i? 2 and R3. The location of the observation 
points is illustrated in Figure In the figures the x-axis is time (in s) and the y-axis is the 
pressure field. 

In the numerical simulation, we assume that the problem is posed in the whole space and to 
solve it, we couple system (|27|l with the perfectly matched absorbing layer model (PML). 
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9.2 Simulation of the response to an incident cylindrical wave of a dike on 
a flexible foundation embedded in a half-space 

To illustrate the efficiency of the method we model in this section the response to an incident 
cylindrical wave of a dike on a flexible foundation embedded in a half-space. This particular 
problem was considered in §32! where it was solved using an expansion of the solution in cylindri- 
cal wave functions in the case of non-dissipative media. In j32| the authors studied this problem 
for different material parameters in order to determine how stiff the foundation should be rel- 
ative to the soil for the rigid foundation assumption in soil-structure interaction models to be 
valid. They concluded that a foundation with the same mass density as the soil but 50 times 
larger shear modulus behaves in rigid manner for this problem. However, for ratios of shear 
moduli less than 16, the rigid foundation assumption is not valid. In this case, soil-structure 
interaction models with a rigid-foundation assumption will not model the differential motion of 
the ground and may underestimate the stresses in the structure (cf. [32J). We consider here a 
ratio of shear moduli equal to 4. Soil-structure interaction is taken into account owing to the 
fact that we discretize the continuous problem. 




free surface 

™ PML 



Figure 13: The geometry of the problem: a dike on a flexible foundation embedded in a half- 
space. 

The geometry of the problem is illustrated in Figure El where T denotes the free surface, f2o 
the hard bedrock, Q\ the flexible foundation and O2 the dike. The physical parameters used in 
the simulation are po = lOOOKg/m 3 , Co = 1450m/s, Qo = +00 in the bedrock, p\ = lOOOKg/m 3 , 
ci = 2900m/s, Qi = 30 in U x and p 2 = 250Kg/m 3 , c 2 = 725m/s, Q 2 = 100 in the dike. The 
angle is equal to tt/2. 

In Figure El we display snapshots of the solution (the pressure field) at different times. 
Diffraction from the free surface is modeled by embedding the solution in a domain of a simple 
shape using the fictitious domain method. To model wave propagation in the infinite half-space 
the fictitious domain is surrounded by an absorbing medium using the PML model. Although 
for this problem a semi-analytical method similar to the one used in |32] can be employed to 
compute the solution, the numerical method proposed in this paper is more general in that it can 
be applied to any complicated geometry and/or propagation media. Moreover, our numerical 
method can be of particular interest in cases where the rigid foundation assumption is not valid 
as it can provide realistic values for the stresses in the structure. 

Conclusion 
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t = 2s 



t = 3s 




Figure 14: Snapshots of the solution: the pressure field in the computational domain at different 
times 
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We employed a rational approximation of the frequency-dependent viscoacoustic modulus in or- 
der to introduce dissipation into time-domain computations. To do so, we followed the approach 
in |18j and chose relaxation frequencies u>/(x) equidistant on a logarithmic scale in the frequency 
range [ w {^ x ; w max ], where w max is the maximal frequency of the used source spectrum. This 
approach will be accurate for propagation in media with a quality factor greater than 10. For 
media with high attenuation (Q < 10) it is necessary in order to obtain accurate results to use 
a non-linear minimization method such as the one proposed in £Q. 

By introducing this approximation of the viscoacoustic modulus into the velocity-pressure 
formulation we obtained a first-order- in-time linear system of equations. To discretize this 
system we used a mixed finite- element method for the discretization in space and a second-order 
finite difference scheme in time. 

The velocity-pressure formulation was coupled with the fictitious domain method in order 
to model the free surface boundary condition on boundaries with complicated geometries, and 
with the PML method to simulate wave propagation in unbounded domains. The efficiency of 
the method was illustrated by numerical results. 
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A Stability analysis 



A.l The continuous problem 

We rewrite the continuous system in time with zero source term, 



<9v 

= Vp, (47) 

%-t a -i-^ («) 

1=1 

dfji 

— +uJir]i = [j, R yidxv{v),\/l. (49) 



By taking the inner products (in I?) of (|4*7]l with v, l|4*%|) with yp — ^^?7; J , and jUjl with rji 

we get 



i=i 



dv \ 

P W V J =(V(p) ' v) ' (50) 

^ ^ _ g^) ' f p_ E^JJ = ^ ^ div ( v )' f^-g^V ( 5i ) 

-^,mj + {^iVhVl) = PRVl (div(v),^) . (52) 



Then, summing (|50|) + — + , we obtain, 



Mjj ~ yi^R 



"I - ') + £ (I (' - 1 ■») • (* - E -») ) + E ^ {% "■) = - E ^ " 



(53) 



Keeping in mind that the energy of the system is, 



1 1 / L L \ L 1 

£ = 2 (/5V ' V) + 2^ [ {P ~ g (P " g J + g 2^ ^' (54) 

we finally get, 

| = -E^IM| 2 <0. (55) 

Which implies that the energy of the system is decreasing with time, when fJ,R and yi are 
positive quantities. The relaxation frequencies u>i are always positive and the same holds for the 
relaxed modulus fiR. The coefficients yi can in practice become negative if we do not solve a 
constraint minimization problem. However, we never encountered in practice a case for which 

-E^IMI 2 >o, 



and thus the problem becomes unstable (in the sense that the energy increases). To avoid this 
instability a constraint minimization algorithm seeking for non-negative yi can be used. 
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A. 2 The discrete problem 

We consider here the more general case where the pressure field is discretized in M\ . Let us 
remark that the discretization space M\ admits the following orthogonal decomposition in L 2 , 

M l h =M h ®{M h ) L , 

where is the space of piecewise constant functions, 

M h = {q h eL 2 /VKe T h , q h \ K € P (K)} , 

and (Mf l ) ± is its orthogonal complement in M\ (with respect to the inner product in I?). To 
simplify the notation, we denote by P the discrete unknown associated with the pressure field 
P = pu € M^, so that we can write P = [Po,P{\ with Po, the projection of P on Mh and Pi the 
projection of P on M^. The memory variables are only discretized on Mh . 
In this case, we can rewrite the discrete system as, (capital letters are used for the discrete 
unknowns and the subscript h is omitted) 

yn+\ _ yn-\ 

P = -B° h P n ~ B\P? (56) 



At 

n+l _ pn L pj"W+l _ fjn 

UrB^V"** (57) 



pn+1 pn L II:. 

±o zJjl _ th zhh. = ll D R°> T v n H 



At ^ At 

i=i 



■n+l _ 

H R B^V n+ i (5* 



pn+l pn 

M ~ M .. n l,T--« ' 1 



At f ~"~' 1 

trn+l rrn w n +l _l / / 

— + w , mviB h 



+ ^ l —= 1 = iJLRyiB^V n+ * (59) 



Then considering the inner products ((|56[) at time (n+l) — f(|56|) at time n)) x V n+ 2 , 

63 x 



get, 



( p o n+1 " + p o™ - £ ' ® x ( P " +1 + P ") ' and ® x ( H i l+1 + ' we 

V i=i i=i J 



p v n+ l : v n+ ^ 



\pV n+ i,V n -^ ) - At [Bl(P% + p n+1 ), y n+ t ) - At [Bl(P? + P? +l ), v n+ ^ 

L 



1=1 

L 



i*> - £ flff + At fx R (B° h ' T V n+1 2 , P£ +1 + iff) - At m B° h ' T V n+1 z , £ H? +1 + fff 



!=1 \ 1=1 

||P™ +1 || 2 = \\P{ 1 \\ 2 + Atfi R {B)f c V n+ ^,P^ +l + P? 

II H n ^^ — H n II 2 

|^" +1 || 2 = ||ifni 2 - w,Ati!-J— g lJ - + MflAtw (5°' T y n +2,^ +1 + Bf 

L 



^ „ . _ © © AB 

Finally summing (JSU) + ^^^^ + + > , we get 



Mfi MR " 2/iA*-R 



£ h ~ £ h _ ST^ \\H l +H l 



At hrvi 



V^- ""' J"" 1 , (64) 
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with the discrete energy being defined by, 

1 i ^ 1 

2e" = fpV^^y^-s) + — ||P™|| 2 + — ||P™ V tf™|| 2 + ||WH| 2 . (65) 



Z=l 



Z=l 



Equation ijMjl shows that the discrete energy is also decreasing, under the same assumptions on 
yi as in IA.1I 

To show under which condition the quantity defined by l|65j) is positive and thus an energy, 
we use the orthogonality relation between Pq and Pi (note that P\ is also orthogonal to Hi), to 
get, 

i L 1 

\H?\\ 2 + 



i 1 

2e n h = ( P {v n+ ^ + v n ~h, (y n+ 5 + v n ~h) + — IIP"!! 2 + V — 



or 



1 



O / ^ \ A+2 



MR 



x M 2 m \\B h f . 2 
4p 



E^r 



where P n = P^ + Pf and P^P™ = P°Po" + B\P™ . We rewrite this equation as a matrix 
associated with the quadratic formulation and we prove that the eigenvalues of this matrix are 
positive under the CFL condition, 



(66) 



with \B%B h \ ^ ^ in ID and \B%B h \ ^ ^ in 2D. 

B Dispersion analysis 

B.l The continuous problem 

Suppose that v(x, t), p(x,t), and t#(x, i)VZ, are plane waves, 

v(x, t) = voexp (i (cut — Kx)), 
p(x, t) = p exp (i (tot - Kx)), 
r/; (x, t) = r/fexp (i (ut - Kx)), 

where Kx = kx in ID and Kx = k x x + k y y = kcos(Q)x + /csin($)y, <E> being the incident angle 
of the plane wave in 2D. Introducing this expression into the time domain system (|23[l. we get 
the dispersion relation, 



U) 



k 2 4 i + E- 



1=1 



1L0 + UJi 



(67) 



with cr = y'p^ ^he relaxed velocity. If the medium is non-dissipative (i.e., y\ = 0V7), l(6"T|l 

becomes the well-known relation uj 2 = K 2 c 2 . Note that the dispersion relation l(6"T|) is no longer 
explicit in lj. 
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B.2 The discrete problem 

We are interested in the general formulation for which the pressure field is discretized in M\ 
and rji in Mh- Considering that V, P, and Hi are plane waves, and employing the same notation 
as in IA.2| we get, 

«rin 2 r-v \ - M2c R In R r, VrO p0 ,r gi^tan ( Xt ) \ . , 

Sm {Xt) ~ "IT [ BhBh + ^ BhBh A^ + 2itan( Xt ) ) (68) 

wherein xt = !± ^, At being the discretization step in time. After some calculations we obtain, 



uAt\ _ A 2 c 2 f^ 2 fk x A x \ 2 (kyAy\\ ( A 2i M taa(aft) \ 



sin 2 I I = — =■ — I sin 



2 J 4 V V 2 ; V 2 J J \ A t ^ + 2itan(^) 

sin 2 f ^ = ^ (sin 2 (**•)) (l + f ^% ) 
V 2 ; 4 V V 2 ^ ^A t ^ + 2itan(^)y 

(69) 

wherein A x and A^ are the discretization step in space. In our case A x = A y = h. 
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